##########################################################################################

library(data.table)
library(optparse)
library(ArchR)

##########################################################################################
option_list <- list(
    make_option(c("--cell_cluster_file"), type = "character"),
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## sperm重新聚类的文件
    cell_cluster_file <- "/public/home/xxf2019/20231121_singleMuti/results/celltype_plot/sperm_recluster/sperm_atac_cluster.tsv"

    ## 整合atac和rna的文件
    comine_data_file <- "/public/home/xxf2019/20231121_singleMuti/results/subcell/cluster2/cluster2.combineRNA.motif_peak2gene.Rdata"

    ## 输出
    out_path <- "/public/home/xxf2019/20231121_singleMuti/results/celltype_plot/sperm_recluster/"

    #comine_data_all_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ_peak-gene/testis_combined_peak.combineRNA.qc.Rdata"


}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_file <- opt$comine_data_file
cell_cluster_file <- opt$cell_cluster_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

a <- load(comine_data_file)
sperm_cell_cluster <- data.frame(fread(cell_cluster_file))

###########################################################################################

c1_cell <- subset(sperm_cell_cluster , cluster == "C1")
c2_cell <- subset(sperm_cell_cluster , cluster == "C2")

###########################################################################################
## ATAC
projHeme5 <- atac_proj
projHeme5@cellColData$sperm_atac_recluster <- "non-sperm"
projHeme5@cellColData$sperm_atac_recluster <- ifelse( rownames(projHeme5@cellColData) %in% c1_cell$cell , "C1" , projHeme5@cellColData$sperm_atac_recluster )
projHeme5@cellColData$sperm_atac_recluster <- ifelse( rownames(projHeme5@cellColData) %in% c2_cell$cell , "C2" , projHeme5@cellColData$sperm_atac_recluster )

corrCutoff <- 0.5       # Default in plotPeak2GeneHeatmap is 0.45

###########################################################################################
## 鉴定C1
sub_proj_c1 <- subsetArchRProject(
  ArchRProj = projHeme5,
  cells = c1_cell$cell,
  outputDirectory = paste0(out_path, "/C1") ,
  dropCells = TRUE,
  force = TRUE
)
projHeme5_tmp <- sub_proj_c1
projHeme5_tmp <- addPeak2GeneLinks(ArchRProj = projHeme5_tmp , useMatrix = "GeneExpressionMatrix")
p2gMat <- plotPeak2GeneHeatmap(
  projHeme5_tmp, 
  corCutOff = corrCutoff, 
  groupBy="sperm_atac_recluster",
  nPlot = 1000000, returnMatrices=TRUE, 
  k=25, seed=1)

kclust_df_out <- cbind( kclust=p2gMat$ATAC$kmeansId , p2gMat$Peak2GeneLinks )
tmp_peak_dat <- data.frame(projHeme5_tmp@peakSet)
tmp_peak_dat$peak <- paste0( tmp_peak_dat$seqnames , ":" , tmp_peak_dat$start , "-" , tmp_peak_dat$end )
kclust_df_out <- merge( kclust_df_out , tmp_peak_dat , by = "peak" )

out_file <- paste0(out_path, "/sperm.reclusters.rna_peakToGene.C1.tsv")
write.table(kclust_df_out , out_file , row.names = F , sep = "\t")

###########################################################################################

sub_proj_c2 <- subsetArchRProject(
  ArchRProj = projHeme5,
  cells = c2_cell$cell,
  outputDirectory = paste0(out_path, "/C2") ,
  dropCells = TRUE,
  force = TRUE
)
projHeme5_tmp <- sub_proj_c2
projHeme5_tmp <- addPeak2GeneLinks(ArchRProj = projHeme5_tmp , useMatrix = "GeneExpressionMatrix")
p2gMat <- plotPeak2GeneHeatmap(
  projHeme5_tmp, 
  corCutOff = corrCutoff, 
  groupBy="sperm_atac_recluster",
  nPlot = 1000000, returnMatrices=TRUE, 
  k=25, seed=1)

kclust_df_out <- cbind( kclust=p2gMat$ATAC$kmeansId , p2gMat$Peak2GeneLinks )
tmp_peak_dat <- data.frame(projHeme5_tmp@peakSet)
tmp_peak_dat$peak <- paste0( tmp_peak_dat$seqnames , ":" , tmp_peak_dat$start , "-" , tmp_peak_dat$end )
kclust_df_out <- merge( kclust_df_out , tmp_peak_dat , by = "peak" )

out_file <- paste0(out_path, "/sperm.reclusters.rna_peakToGene.C2.tsv")
write.table(kclust_df_out , out_file , row.names = F , sep = "\t")